---
title: "Tityus achillis venom spraying 2024"
author: "Léo Laborieux"
date: "2024-04-23"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=TRUE}
library(ggplot2)
library(dplyr)
library(car)
library(lme4)
library(coin)
library(emmeans)
library(brms)
library(bayestestR)
library(see)
```

## Preprocessing

Reading data, setting classes

```{r Preprocessing}
df <- read.csv("/Users/leo/Documents/MEME/Tityus Project/Final script/tityusdata_LL2024.csv", sep=",")
df <- as.data.frame(apply(df, 2, function(x) ifelse(x == "", NA, x)))
df$Type <- as.factor(df$Type)
df$Velocity <- as.numeric(df$Velocity)
df$Initial.angle <- as.numeric(df$Initial.angle)
df$Initial.height <- as.numeric(df$Initial.height)
df$Angle.min <- as.numeric(df$Angle.min)
df$Angle.max <- as.numeric(df$Angle.max)
df$Specimen <- as.integer(df$Specimen)
df$Duration.frames <- as.integer(df$Duration.frames)
df$Specimen <- as.factor(df$Specimen)
df$Duration.frames[df$Type=="Flick"] <- 1

```

## Individual data

```{r Venom ducts + ind data}


ind_data <- read.csv("/Users/leo/Documents/MEME/Tityus Project/Final script/ind_data.csv", sep=";")
to_estimate <- data.frame(Prosoma.L = ind_data$P.length)
ind_data$Number <- as.integer(ind_data$Number)

matching_rows <- match(df$Specimen, ind_data$Number)

# Add Sex, Prosoma length and Duct values to df based on matching rows
df$Sex <- ind_data$Sex[matching_rows]
df$Sex <- as.factor(df$Sex)
df$Prosomalength <- ind_data$P.length[matching_rows]

```

## Calculation of ranges (equation 2)

```{r Ranges (eq2)}

df$Range <- df$Velocity*cos(df$Initial.angle*(pi/180))*(1/9.81)*(df$Velocity*sin(df$Initial.angle*(pi/180)) + sqrt((df$Velocity*sin(df$Initial.angle*(pi/180)))^2+2*9.81*(df$Initial.height/100)))

```

## Calculation of apexes (equation 3)

```{r Apexes (eq4)}

df$Apex <- df$Initial.height/100 + (df$Velocity)^2 * (sin(df$Initial.angle*(pi/180)))^2/(2*9.81)

#in meters
```

## Calculation of secreted volumes (equation 4)

```{r Secreted volumes (eq1)}

df$Volume <- pi*(1.25*10^-4)^2* df$Velocity * (df$Duration.frames / 240)
df$Volume <- df$Volume*1000

#in L
```

## Plots

```{r Plots, echo=FALSE}

#####theme ------

require(grid)
pub<- theme_update(
  panel.grid.major=element_line(colour=NA),
  panel.grid.minor=element_line(colour=NA),
  panel.background = element_rect(colour = NA,fill=NA,size=1),
  panel.border = element_rect(linetype = "solid", colour = "gray47",fill=NA),
  axis.line.x = element_line(color="black"),
  axis.line.y = element_line(color="black"),
  axis.title.x=element_text(size=12,face="bold",hjust=0.5,vjust=0.5,angle=0),
  axis.title.y=element_text(size=12,face="bold",hjust=0.5,vjust=0.5,angle=90),
  axis.text.x=element_text(colour="black",angle=0,size=12),
  axis.text.y=element_text(colour="black",angle=0,size=12),
  axis.ticks=element_line(colour="black",size=1))

##### boxplot velocities -----
ggplot(df, aes(x = Type, y = Velocity, fill = Type)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.6, aes(size=Duration.frames)) +
  scale_size(range = c(2, 10)) +
  scale_fill_manual(values = c("ivory4", "salmon2")) + 
  labs(x = "Pulse type", y = "Velocity") +
  theme(legend.position = "bottom",
        plot.margin = margin(30,30,30,30))


##### boxplot hi/L -----

ggplot(df, aes(x = Type, y = Initial.height/(Prosomalength/10), fill = Type)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.6) +
  scale_size(range = c(2, 10)) +
  scale_fill_manual(values = c("ivory4", "salmon2")) +
  labs(x = "Pulse type", y = "hi/L") +
  theme(legend.position = "bottom",
        plot.margin = margin(30,30,30,30))

##### boxplot range -----

ggplot(df, aes(x = Type, y = Range, fill = Type)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.6) +
  scale_size(range = c(2, 10)) +
  scale_fill_manual(values = c("ivory4", "salmon2")) +
  labs(x = "Pulse type", y = "Range (m)") +
  theme(legend.position = "bottom",
        plot.margin = margin(30,30,30,30))


##### angular distributions ------
ggplot(df, aes(x = Initial.angle, fill = Type, color = Type)) +
  geom_density(alpha = 0.7) +
  labs(x = "Angle", y = "Density", title = " ") +
  xlim(-180,180) +
  coord_polar(start=pi/2) + 
  theme_minimal() +
  scale_fill_manual(values = c("ivory4", "salmon2")) +
  scale_color_manual(values = c("ivory4", "salmon2")) +
  theme(panel.grid.minor = element_blank(),
        axis.text.x = element_text(margin = margin(t = 10)),
        axis.text.y = element_text(margin = margin(r = 10)),
        legend.position = "none",  # Remove legend for "Type"
        plot.margin = margin(30,30,30,30))




##### boxplot volume ------
ggplot(df, aes(x = Type, y = log10(Volume), fill = Type)) + #### log transformation for readability
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.08, alpha = 0.6) +
  scale_size(range = c(2, 10)) +
  scale_fill_manual(values = c("ivory4", "salmon2")) +
  labs(x = "Attack type", y = "Volume expelled (L)") +
  theme(legend.position = "bottom",
        plot.margin = margin(30,30,30,30))



##### trajectories (equation 2) ------

  
# Define the ballistic equations in terms of x and y
ballistic_equation <- function(x, angle, velocity, y0) {
  gravity <- 9.81  # Acceleration due to gravity in m/s^2
  y <- y0 + tan(angle) * x - gravity * (x^2/(2*velocity^2*cos(angle)^2))
  return(data.frame(x = x, y = y))
}

# Generate x values for plotting
x_values <- seq(0, 0.4, by = 0.001) 

# Create a df to store trajectory data
all_trajectories <- data.frame(x = numeric(0), y = numeric(0), type = factor(character(0)), color = character(0))

# Plot specific trajectories from dataframe df
for (i in 1:nrow(df)) {
  y0 <- df$Initial.height[i] / 100  # Convert to meters
  velocity <- df$Velocity[i]
  angle <- df$Initial.angle[i] * pi / 180  # Convert to radians
  type <- df$Type[i]
  
  # Compute y values for the current trajectory
  trajectory_data <- ballistic_equation(x_values, angle, velocity, y0)
  
  # Add trajectory data to the dataframe
  trajectory_data$type <- factor(paste("Trajectory", i))
  trajectory_data$color <- type  # Assign color based on df$Type
  
  all_trajectories <- rbind(all_trajectories, trajectory_data)
}

# Plot all trajectories
ggplot(all_trajectories, aes(x = x, y = y, group=as.factor(type), color = color)) +
  geom_line() +
  xlim(-0.05, 0.4) +
  ylim(0, max(all_trajectories$y)) +
  labs(x = "Horizontal Distance (m)", y = "Vertical Distance (m)", title = "Ballistic Trajectories") +
  scale_color_manual(values = c("ivory4", "salmon2")) +  # Set colors for the two levels of df$Type
  coord_fixed() +
  theme_minimal()

```

## Statistical analyses

```{r Stats}

##### initial velocities ----

#linear model

lm_velocity <- lmer(Velocity ~ Type + (1|Specimen), data=df) #singular
lm_velocity <- lm(Velocity ~ Type, data=df)
summary(lm_velocity)
Anova(lm_velocity)

emmeans(
  object = lm_velocity,
  specs = pairwise ~ Type,
  regrid = "response"
)

##### hi/L ----


#linear model

lm_hi <- lmer(Initial.height ~ Type + (1|Specimen), data=df) #not singular
summary(lm_hi)
Anova(lm_hi)

emmeans(
  object = lm_hi,
  specs = pairwise ~ Type,
  regrid = "response"
)

##### ranges ----

#linear model

lm_range <- lmer(Range ~ Type + (1|Specimen), data=df) #singular
lm_range <- lm(Range ~ Type, data=df)
summary(lm_range)
Anova(lm_range)

emmeans(
  object = lm_range,
  specs = pairwise ~ Type,
  regrid = "response"
)


##### apexes ----

#linear model

lm_apex <- lmer(Apex ~ Type + (1|Specimen), data=df) #singular
lm_range <- lm(Apex ~ Type, data=df)
summary(lm_apex)
Anova(lm_apex)

emmeans(
  object = lm_apex,
  specs = pairwise ~ Type,
  regrid = "response"
)

#prcc

##### initial angles ------

# Mixed-effects model with circular distribution using brms
fit <- brm(
  bf((Initial.angle*3.141593/180) ~ Type  + (1 | Specimen)),  #model specification after converting to radians
  data = df,
  family = von_mises(),  # Specify the circular distribution
  prior = set_prior("normal(0,(3.141593/6))", class = "b"),  # Specify a normal (mu=0, sd=pi/6) prior
  iter=42000,
  warmup=2000
)

summary(fit)
plot(fit, variable = c("b_TypeSpray"))

# 
# test_nullregion <- bayesfactor_parameters(fit, null = c(-0.1, 0.1))
# test_nullregion
# plot(test_nullregion)


# bayes-factor with directional hypothesis ("sprays are initiated at a higher angle than flicks")
test_directional <- bayesfactor_parameters(fit, direction = ">")
test_directional
plot(test_directional)


#####volumes ----

#linear model
lm_vol <- lmer(1000000*Volume ~ Type + (1|Specimen), data=df) #not singular
summary(lm_vol)
Anova(lm_vol) 

emmeans(
  object = lm_vol,
  specs = pairwise ~ Type,
  regrid = "response"
)

```

## Misc calculations

```{r Other values reported in text}

#####calculating the maximum volume expelled by a single specimen during a single session -----
max_expelled_volume <- sum(df$Volume[df$Round == 1 & df$Specimen == 5])
print(paste("Maximum volume expelled by specimen 5 during Round 1:", max_expelled_volume))

#####general info on spray durations ------
spray_duration_summary <- summary(df$Duration.frames[df$Type=="Spray"])
print("Summary of spray durations for 'Spray':")
print(spray_duration_summary)

spray_duration_sd <- sd(df$Duration.frames[df$Type=="Spray"])
print(paste("Standard deviation of spray durations for 'Spray':", spray_duration_sd))

#####means and standard deviations for velocity, range, apex, and volume per pulse type ------
#velocity
mean_velocity_spray <- mean(df$Velocity[df$Type=="Spray"])
sd_velocity_spray <- sd(df$Velocity[df$Type=="Spray"])
print(paste("Mean velocity for 'Spray':", mean_velocity_spray))
print(paste("Standard deviation of velocity for 'Spray':", sd_velocity_spray))

mean_velocity_flick <- mean(df$Velocity[df$Type=="Flick"])
sd_velocity_flick <- sd(df$Velocity[df$Type=="Flick"])
print(paste("Mean velocity for 'Flick':", mean_velocity_flick))
print(paste("Standard deviation of velocity for 'Flick':", sd_velocity_flick))

#range
mean_range_flick <- mean(df$Range[df$Type=="Flick"])
sd_range_flick <- sd(df$Range[df$Type=="Flick"])
print(paste("Mean range for 'Flick':", mean_range_flick))
print(paste("Standard deviation of range for 'Flick':", sd_range_flick))

mean_range_spray <- mean(df$Range[df$Type=="Spray"])
sd_range_spray <- sd(df$Range[df$Type=="Spray"])
print(paste("Mean range for 'Spray':", mean_range_spray))
print(paste("Standard deviation of range for 'Spray':", sd_range_spray))

#apex
mean_apex_flick <- mean(df$Apex[df$Type=="Flick"])
sd_apex_flick <- sd(df$Apex[df$Type=="Flick"])
print(paste("Mean apex for 'Flick':", mean_apex_flick))
print(paste("Standard deviation of apex for 'Flick':", sd_apex_flick))

mean_apex_spray <- mean(df$Apex[df$Type=="Spray"])
sd_apex_spray <- sd(df$Apex[df$Type=="Spray"])
print(paste("Mean apex for 'Spray':", mean_apex_spray))
print(paste("Standard deviation of apex for 'Spray':", sd_apex_spray))

#volume
mean_volume_flick <- mean(df$Volume[df$Type=="Flick"])
sd_volume_flick <- sd(df$Volume[df$Type=="Flick"])
print(paste("Mean volume for 'Flick':", mean_volume_flick))
print(paste("Standard deviation of volume for 'Flick':", sd_volume_flick))

mean_volume_spray <- mean(df$Volume[df$Type=="Spray"])
sd_volume_spray <- sd(df$Volume[df$Type=="Spray"])
print(paste("Mean volume for 'Spray':", mean_volume_spray))
print(paste("Standard deviation of volume for 'Spray':", sd_volume_spray))

#####mean and standard deviation of angular range during sprays
mean_angular_range <- mean(abs(df$Angle.min - df$Angle.max), na.rm = TRUE)
sd_angular_range <- sd(abs(df$Angle.min - df$Angle.max), na.rm = TRUE)
print(paste("Mean angular range during sprays:", mean_angular_range))
print(paste("Standard deviation of angular range during sprays:", sd_angular_range))

##### volume ratio between the largest spray and the smallest 
difference_vol <- max(df$Volume) / min(df$Volume)
print(paste("Volume ratio between the largest and smallest spray:", difference_vol))

##### estimation of specimen's 7 "prevenom" reserve
vol_before_white <- sum(df$Volume[df$Round==3 & df$Specimen==7])
print(paste("Estimation of specimen 7's 'prevenom' reserve during Round 3:", vol_before_white))

```
